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Auto and crosscorrelograms for the spike response of LIF neurons with slow synapses 
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An analytical description of the response properties of simple but realistic neuron models in the 
presence of noise is still lacking. We determine completely up to the second order the firing statistics 
of a single and a pair of leaky integrate-and-fire neurons (LIFs) receiving some common slowly filtered 
white noise. In particular, the auto- and cross-correlation functions of the output spike trains of 
pairs of cells are obtained from an improvement of the adiabatic approximation introduced in [lj. 
These two functions define the firing variability and firing synchronization between neurons, and 
are of much importance for understanding neuron communication. 
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The variability of the spike trains of cortical neurons 
and their correlations might constraint the coding ca- 
pabilities of the brain 0, but they can also reflect the 
strategies the brain uses to decipher the stimuli arriving 
from the world Q . Neurons in cortex fire with high vari- 
ability resembling Poisson spike trains and nearby 
pairs of cortical neurons fire in a correlated fashion [2j, 
reflecting the presence of some common source of noise. 
These variability and correlation of the spike trains af- 
fect the firing statistics of a neuron receiving those inputs 
It has been shown that the large variability ob- 
served in vivo can be accounted for by neuron models 
operating in a regime in which the membrane time con- 
stant, r TO , becomes shorter or comparable to the synaptic 
decay constants, T s ,due to spontaneous background ac- 
tivity (t s > T m ) However, very little progress has 
been made in providing analytical tools to describe such 
variability and correlations found in cortex. 

In this Letter we study analytically the variability and 
correlations in the firing responses of pairs of LIF neurons 
receiving both common and independent sources of white 
noise input filtered by synapses in the regime t s > r m . 
For a single neuron we obtain the firing rate, the auto- 
correlation function of its output spike train (ACF), the 
Fano factor of the spike count, F^. For a pair of cells, we 
obtain the crosscorrelation function of their output spike 
trains (CCF) and the correlation coefficient of their spike 
counts, p. These results characterize completely the fir- 
ing response of these spiking neurons up to second order, 
and open the possibility for a principled way of including 
synchrony effects in the modeling of biologically plausible 
spiking neural networks. 

The neuron and input models. The membrane po- 
tential V(t) of a single LIF neuron with membrane time 
constant r m and receiving an afferent current I(t) obeys 



9] which is filtered by synapses with decay time constant 
Ts, resulting in a current described by 



r s i{t) = -/(*) + f/, + ar}{t) , 



(2) 



where r)(t) is a Gaussian white noise with zero mean and 
unit variance. We simplify eqs. by performing the 

linear transformations I = fi + z <j/\/2t s and V = /j,T m + 
x o~^J T m /2, obtaining 
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with 7 = y/T m /T s . In the normalized potential, x, the 
threshold and reset read 8 = \/2(0 — /UT TO ) /o~-JT m and 

H = \/2(H - flT m )/(Ty/n~. 

The autocorrelation function. To determine the 
ACF, first we describe the time evolution of the prob- 
ability density of having the neuron in the state (x, z) 
at time t given that initially the neuron has just fired 
{x = H) and z — zq. The Fokker-Planck equation (FPE) 
for this density, P(x, z, t\H, zq), is [l| 
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where e = 7= y/T m /r s and L z = -§^z + -§^- J(z,t\z ) is 
the probability density of having a spike at time t along 
with a fluctuation z given that z — zq at time t = 0. This 
probability is expressed as a function of the density P as 
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A spike is generated when V{t) reaches a threshold 8, 
after which the neuron is reset to H, from where it con- 
tinues integrating the current @. The external input is 
modeled by a white noise with mean \jl and variance a 2 



J(z,t\z ) 



-(-Q + -yz)P(e,z,t\H,zo). (6) 



Solving the FPE |5]| with J(z,t\z ) as a source term at 
x = H means that each time a spike is produced, the 
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FIG. 1: ACF of a LIF neuron in the sub (left) and 
suprathreshold regimes (right). The figures show the typi- 
cal shape in both regimes: no oscillations and a large peak 
in the subthreshold regime (fiT m < 0) and damped oscilla- 
tions in the suprathreshold regime (/ir m > O). Thick lines 
are the analytical results obtained from eq. (|12[l (the sum has 
been cut at re = 200 with t = 200ms) and thin lines corre- 
spond to the numerical simulations of the same LIF neuron. 
Parameters for the subthreshold (suprathreshold)neuron are 
H = 85Hz(U5Hz) a 2 — 6Hz(3Hz). Other parameters are 
H = 0, G = 1, T m = 10ms and t 3 = 20ms. 



normalized potential x is reset to H while z keeps its 
same value. 

The integral J dzJ(z,t\zo) expresses the probability of 
having a spike at time t conditioned to the fact that 
z — z at time t = 0. We define the ACF, C(t), as the 
probability density of firing a spike at time t > condi- 
tioned to the fact that at time t = there was a spike. 
Therefore, C(t) is the average of J dzJ(z,t\zo) with the 
distribution of Zq conditioned to the production of a spike 
at time t = 0, B(zq). Since B(z) is the distribution of 
z at the moment of a spike, then B(z) = J(z)/v, where 
J(z) is the limit t — > oo of J(z, t\zo), and v is its normal- 
izing factor (y = j dzJ(z)) and also the firing rate of the 
LIF neuron defined by eqs. ((3j{4j) . Therefore, the ACF is 
computed as 
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The solution of the FPE ^ and eq. is simplified by 
noticing that z is a pure Ornstein-Uhlenbeck process, eq. 

and therefore its marginal distribution, P(z,t\zo), is 
(see, e.g., [§]) 
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which broadens over time and for ( > r s approaches a 
normal distribution, p(z) = e~ z / 2 /\/2tt. 
The analytical solution. We expand both 
P(x, z, t\H , Zq) and J(z,t\z Q ) in powers of e 2 , as P = 
P a + e 2 Pi + 0(e 4 ) and J = J + e 2 J l + 0(e 4 ), following 



a technique introduced in [l[ for the stationary FPE. In 
this expansion, the parameter 7 in eqs (|5I6[) is assumed 
to be fixed. Only at the end, when the leading orders of 
the expansion have been found, 7 is given its true value 
7 = \Ji~ m lT s . 

The solution at zero-th order of the FPE (O satisfying 
conditions (|6l8p is 

P (x, z, t\H, zq) - P(z, t\z )6(x - X(z, t)) , (9) 

where X(z, t) is the time evolution of the variable x ob- 
tained from eq. |3]) with frozen z and initial condition 
H. Notice that x = X(z,t) is a periodic function of t, 
because whenever x — 0, x is reset to H. Its period, 
T(z) = r m ln(H - 7 z/e - 7 z) (T(z) = 00 for z < 6/7), 
is the inter-spike-interval (ISI) of a LIF neuron receiving 
a frozen z, and it is calculated from eq. ([3]) as the first 
time T at which X(z, T) — 0. After expressing the delta 
functions in terms of t, the probability density current, 
eq. ([6]), at zero-th order becomes 



J {z, t\z ) = P{z, t\z ) S(t - nT{z)) . 



(10) 



This expression has a simple interpretation. The sum of 
delta functions in the index n represents a regular train 
of spikes with ISI T(z), as if z were fixed. Therefore, the 
probability of having a spike along with a fluctuation z 
at time t, Jo, is given at a first approximation by the 
product of both the probability of finding at time t a 
spike of the train generated with frozen fluctuation z, 
and the probability of having such a fluctuation z at time 
t starting from the initial condition z = zq, P(z,t\zo). 
Note that in eq. (fTQ| the noise is allowed to evolve in 
time following the distribution P(z, t\zo). It has been 
proved that the stationary (frozen) distribution of z can 
be employed to describe the firing rate of LIF neurons [l|, 
7] , and used the approximation that z is constant during 
the ISIs to describe the Fano factor of non-LIF neurons 
with weak noise However, freezing completely the 

noise z in eq. (|10p leads to very poor predictions in our 
problem (not shown). 

To determine the ACF, eq. (O, at zero-th order, we 
have to the zero-th order J(z) is required, which is [l| 



J (z) = v (z)p{z) , 



(11) 



where vq{z) = 1/T(z) for z > 0/7 and vq{z) — other- 
wise. Finally, the zero-th order ACF is computed, after 
using eqs. (|7I10I11[) and evaluating the delta functions, 
as 



Co{t) = £ 



n=l 



dz Jq{z ) (-jz n - H)(jz n - 0) 
u T m n 7 (0 - H) 



P(z n ,t\z ) 
(12) 
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FIG. 2: The firing rate (left), the Fano factor (right), Fn, and 
the CV 2 (right) for the output spike train of a LIF neuron are 
plotted as a function of r s . The firing rate prediction (line) is 
calculated using both eqs. {13} and Vq = J dzJo(z), and it is 
compared with simulation results (points) . The Fn predicted 
by eqs. (| 14112113ft (line) is compared with the Fn (squares) 
and CV 2 (triangles) obtained from simulations. Parameters 
values are as in the subthreshold regime of Fig. |T} ■ 

where z n = z n (t) = 7- x (6 - ffe~ t /' iT '")/(l - e -*/" T ™). 
The z n s are the roots of the equations t = nT(z n ), the 
zeros of the delta functions in eq. (jTUJ) . 

In Fig. {1} we plot the ACF for the output spike train 
of a LIF neuron computed using eq. (|12p and compare it 
with simulation results. The agreement is very good in 
both the subthreshold (left) and suprathreshold (right) 
regimes. In both regimes, the ACF shows a prominent 
peak after a relative refractory period of about 10ms 
(~ T m ). This means that the potential has to be in- 
tegrated from reset to threshold to emit the first spike. 
The prominent peak indicates that the neuron is bursty, 
producing spikes that are grouped within short time in- 
tervals of 20ms (~ r s ) [![. After the prominent peak, 
the ACF decays to a steady-state value either monotoni- 
cally (left) or with a damped oscillation (right). Damped 
oscillations are a robust feature in the suprathreshold 
regime, as is their absence in the subthreshold regime. 
This reflects the fact that the neuron in the suprathresh- 
old regime fires more regularly, and therefore the output 
spikes tend to occur at integer number of times the mean 
ISI (see the peaks of the oscillations in the ACF). For long 
times (t > t s ) the memory of the spike at time t = has 
disappeared, and the ACF decays to the unconditioned 
probability of having a spike, that is, the firing rate of 
the LIF neuron. 

The firing rate, Fano factor and CV. As it is clear, 
the firing rate can be obtained from the ACF, eq. (fT2")> . 
in the limit of long times (t ^S> t s ). This rate has the 
expression 

v = lim V hznjt) - H)(jz n (t) - Q) e - z „ {t f/2 
° sphtTm n 7(6 — H) 

(13) 

A different expression for the firing rate can be computed 



FIG. 3: Left: Theoretical (thick line, eq. (|17[0 and simulated 
(thin line) CCFs normalized by the firing rate of one of the 
neurons as a function of time lag. Here a 2 — 2Hz. Right: 
Theoretical (line, eq. (|18[l 1 and simulated (points) correlation 
coefficient, p, for the output spike trains of a pair of identical 
LIF neurons as a function of the fraction of common noise, 
a 2 la 2 . The numerical p is calculated using eq. (|18[1 inte- 
grating the simulated CCF. Parameters for both figures are 
p = 85Hz, a 2 — 9Hz, and the others as in Fig. ([lj. 

using vq = J dzjg(z) [l[. In fact, both expressions give 
identical results when they are plotted as a function of r s 
(continuous curve in Fig. ([5]), left). However, computa- 
tionally, eq. (|13[) is much faster because it only involves a 
sum that can be cut at n ~ 200 (using t — 200ms). Nat- 
urally, the number of terms needed to approximate the 
ACF and the firing rate grows as t increases. Comparison 
of both expressions of vq with simulation results shows 
that the prediction is very good even when t s ~ r m . 

The Fano factor of the output spike train, Fn, defined 
as the ratio between the variance of the spike count and 
its mean evaluated for long time windows, is directly re- 
lated to the time integral of ACF as ([ll| and see, e.g., 
eq. (3) of ref @) 

POO 

F N = 1 + 2 dt (C(t) - v) . (14) 
Jo 

We have evaluated the zero-th order Fjy in eq. (fL4"| us- 
ing the zero-th order solutions of C(t) and v, eqs. (|12U3p . 
The prediction fits very well the simulation results (right 
panel of Fig. (|2|)). We have also computed the coeffi- 
cient of variation of the ISIs, CV, of the neuron response 
using simulations (same panel). It is known that for re- 
newal processes F N = CV 2 (e.g. for a Poisson process 
F N ee CV 2 = 1, and C(t) = v). Here we find that 
Fn ~ CV 2 even when the output response is not a re- 
newal process. This is because, although the synaptic 
time scale introduces correlations in the successive ISIs, 
since for low (but typical) rates r s < v~ l , the correla- 
tion between successive ISIs is small. This explains the 
similarity between Fn and CV 2 . Notice that the firing 
variability is large when r s > r m 0, HJ . 
The crosscorrelation function and correlation co- 
efficient. A central issue to describe population dynam- 
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ics is to understand the way neuron activity synchronizes. 
Here we study a pair of identical LIF neurons (k = 1,2) 



-V k +T m ( I k (t) + I c {t) ) 



(15) 



receiving both an independent source of current, I k (t), 
and a common source, I c {t). Each current is described 
by an equation identical to eq. (JSJ), with mean fj,i n d and 
variance o 2 nd for the independent components, and mean 
/i c and variance a 2 for the common component. Each 
neuron receives a total mean current [i = /Uj ra( j + and 
total variance a 2 = a 2 nd + a 2 . 

The CCF of the output spike trains of the two neu- 
rons (denoted as CC(A)) can be obtained by an anal- 
ysis similar to that used for the ACF. The CCF is de- 
fined as the joint probability density of having a spike 
of neuron 1 at a given time and a spike from neuron 
2 after a delay A. Here we only summarize the main 
results. First, we define the normalized fluctuations 
u k = (Jjt+7 C — /i)/<7, having zero mean and unit variance. 
Notice that these are not independent because of the 
common input I c . Second, if neuron 1 has a fluctuation 
ui, the probability density that after a delay A neuron 
2 has a fluctuation u 2 , P(u2, A\ui), is a Gaussian dis- 
tribution with mean (u2(A,u\)) = u\e^ A l Ta a 2 ja 2 and 



variance Var(u2(A)) 
long t s 



1 - e- 2A / T °a 2 c /o- 2 . Then, for 



where a n = z n (t) and b m = z rn (t + A), with z n (t) as 
in eq. ([12")) . The theoretical CCF matches very well the 
simulated one (Fig. J3]), left). Typically, the prediction 
underestimates the central peak occurring at time lag 
zero. The central peak of the CCF decays within a time 
of the order of t s (notice that the CCF is symmetric 
around A = 0). This is because the synaptic input, being 
slower than the neuron dynamics, sets its own time scale 
in the dynamics of interactions of the two neurons. The 
existence of a single peak is robust for low values of a 2 in 
both the sub and suprathreshold regimes, but other side 
secondary peaks arise when all the noise is essentially 
common. For long A, the CCF converges to the product 
of the firing rates at zero-th order, v 2 (see eq. (US])), 
because the neurons fire independently. 

The correlation coefficient, p, of the spike counts for 
long time windows of the output spike trains of two iden- 
tical neurons can be computed from their CCF ([ll[ and 
see, e.g., eq. (4) of ref [(|) 



F N v 



ds {CC(s) - v 2 ) . 



(18) 



CC (A) = lim / duidu 2 P(u 2 , A\ui) p(m) 

t->oaJ 

oo 

5{t-nT{ux)) S(t + A-mT(u 2 )) , (16) 

n,m— 1 

where T(u k ) = r m ln(H - ju k /Q- ju k ) {T(u k ) = oo for 
u k < /j) is the ISI of the neuron i receiving a constant 
fluctuation u k , and p{u\) is a normal distribution de- 
scribing the steady state distribution of the fluctuations 
of neuron 1. The quantities 7, H and are defined as 
before. The two sums of delta functions in eq. p^|) can 
be interpreted as the product of two output spike trains 
with fixed ISI (determined by the input fluctuations), 
quantity which has to be averaged over all the possible 
fluctuations. The result of such an average is the CCF 
when the limit t — > 00 is taken to allow randomization 
of the initial conditions, cq. p6p . This equation can be 
simplified by integration of the delta functions, obtaining 



CC (A) = lim V 



(7C - H)( 7 b n - H) 



^ n m t^j 7 2 (0 — H) 2 
{ia m - 0)( 7 6 m - 0)P(6 m , A\a n ) p(a n ) , (17) 



For the two neurons in eq. (|15p it can be computed at 
zero-th order using the zero-th orders of CC(A), eq. (fT"7|) . 
Fn and v. We have compared the theoretical and sim- 
ulated p as the fraction of common noise increases (Fig. 
©, right). The prediction is good for low values of com- 
mon noise, and departs from the simulations for larger 
values. As the common noise increases, p increases mono- 
tonically and reaches p — 1 when the common noise 
equals the total input noise. Correlation coefficients of 
~ 0.1 as those found in cortex are predicted accu- 
rately, and they are obtained when the common noise 
represents ~ 20 per cent of the total synaptic noise en- 
tering into the neuron, which can be a realistic value (2j. 
Therefore, the right plot at Fig. 3 provides a valuable tool 
to estimate the fraction of common noise from the corre- 
lations of the spike trains of pairs of neurons, a quantity 
which otherwise is not available experimentally. 

The results we have obtained at the cell level open the 
way for a systematic investigation of the role of correla- 
tions in neuronal networks. 
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